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We study the dynamical instabilities of superfluid flows in the 5=1 Bose-Hubbard model. The time evo¬ 
lution of each spin component in a condensate is calculated based on the dynamical Gutzwiller approximation 
for a wide range of interactions, from a weakly correlated regime to a strongly correlated regime near the Mott- 
insulator transition. Owing to the spin-dependent interactions, the superfluid flow of the spin-1 condensate 
decays at a different critical momentum from a spinless case when the interaction strength is the same. We 
furthermore calculate the dynamical phase diagram of this model and clarify that the obtained phase boundary 
has very different features depending on whether the average number of particles per site is even or odd. Finally, 
we analyze the density and spin modulations that appear in association with the dynamical instability. We find 
that spin modulations are highly sensitive to the presence of a uniform magnetic field. 
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I. INTRODUCTION 

A recent experimental development in optical lattices offers 
the unprecedented potential to study the dynamical properties 
of many-body interacting ultracold atoms mu. In particu¬ 
lar, the superfluid flow of a Bose-Einstein condensate (BEC) 
loaded on a lattice exhibits a novel instability called dynam¬ 
ical instability that was predicted a decade ago |3, kj], which 
has attracted much attention both theoretically |5 <S| and ex¬ 
perimentally mu. The dynamical instability is induced by 
the interplay between the lattice periodicity and nonlinearity 
due to the inter-particle interactions in the BEC. When the 
system becomes dynamically unstable, the energy of an ex¬ 
cited mode has an imaginary part 0 ). Therefore, an arbitrary 
small density fluctuation in a uniform superfluid flow grows 
exponentially in time, resulting in a drastic decay of the orig¬ 
inal flow. These features are in contrast with the well-known 
Landau instability, which is the energetic instability caused by 
decaying from the initial metastable state. 

Dynamical instability itself is widely seen in various non¬ 
linear systems governed by classical fluid mechanics. How¬ 
ever, using ultracold atoms, we can now advance the study of 
superfluid instabilities to the next stage, namely dynamical in¬ 
stabilities in systems with internal degrees of freedom. It has 
been known that multicomponent systems exhibit rich physics 
such as diverse quantum phases in an optical lattice lfl4l420h . 
and the dynamical instability of multicomponent bosons has 
also recently been studied KlHMI . Moreover, bosons with 
unfrozen spin degrees of freedom specifically exhibit com¬ 
plex and intriguing phenomena caused by spin mixing pro¬ 
cesses tf25l l. The spin-1 bosons have therefore been inves¬ 
tigated intensively as the simplest bosonic system with un¬ 
frozen spin degrees of freedom. A series of studies have 
revealed interesting instabilities in the spin-1 BEC based on 
the Gross-Pitaevskii equation, spin mixing instability lf26ll27tl . 
spin counterflow instabil ity I28ll29il . and the spontaneous for¬ 
mation of spin domains |j3( iw- These phenomena are spe¬ 


cific to the spin-1 bosonic system and have hardly been under¬ 
stood only by conventional linear stability analysis GHIH 
because of the spin mixing process inherent in the system. 

The spin-1 bosons in optical lattices are well described by 
the 5 = 1 Bose-Hubbard model (BHM) ll32l - l34}l . The phase 
diagram and the static properties of this model have already 
been extensively studied using several theoretical methods 
ll35l - l43ll . From these studies, the Gutzwiller-type variational 
wave functions are good at capturing the superfluid (SF) to 
Mott-insulator (MI) transition in the 5 = 1 BHM, aside from 
spin correlations in the MI phase ll43l -(45ll. It has been found 
that the SF-MI transition in this model strongly depends on 
whether the average particle number per site, n, is even or 

odd mm. 

We study the effect of spin interaction in the dynamical 
instability according to the following three interests. First, 
whether does the parity about the average number of parti¬ 
cles per site as mentioned above also appear in the dynamical 
phase diagram or not? This motivates us to explore a role 
of spin degrees of freedom in dynamical phenomena of a su¬ 
perfluid, which remains to be clarified. Second, how do spin 
mixing processes among spin components in spin-1 superfluid 
flows affect the dynamical instability? Spin mixing, which 
is an important feature in bosonic spin systems, does not ex¬ 
ist in classical fluids and multicomponent cold atom systems. 
Therefore, the effect of spin mixing on instabilities of fluids 
itself is intriguing. Finally, we are interested in the very re¬ 
cent development of experimental techniques for observing 
spin dynamics of condensates in optical lattices as reported 
by L. Zhao et al. [46j|. It was revealed experimentally that the 
intensity of lattice potential significantly affects spin mixing 
dynamics in a spin-1 system. The experiment of the dynam¬ 
ical instability in a spin-1 system is therefore expected to be 
demonstrated in the near future. 

In this paper, we analyze the dynamical instability of the 
spin-1 condensate in the 5 = 1 BHM for a wide range of in¬ 
teraction parameters with antiferromagnetic or ferromagnetic 



2 


interactions, focusing on the stability of spin-resolved super¬ 
fluid flows. First, we reveal how the spin mixing process af¬ 
fects the real-time evolution of each spin component in the 
flow. We employ the dynamical Gutzwiller approximation 
that was used by Altman et al. 17. 81 to analyze the dynamical 
instability of an SF in the spinless BHM. Recently, Natu et 
al. also applied this method to the 5=1 BHM and they cal¬ 
culated the low-lying excitation spectrum [||47]. We show the 
dynamical decay of the 5=1 superfluid flow and the corre¬ 
sponding time development of the spin components. Next, we 
demonstrate the parity dependence of dynamical instability in 
the 5=1 BHM constructing dynamical phase diagrams. In 
the antiferromagnetic case, the stable flow region on the phase 
diagram shrinks when the average number of particles is odd, 
while it grows for the even average numbers compared with 
the spinless case. We find that this phenomenon is caused 
by the spin mixing process. Finally, we discuss the density 
and spin modulations associated with the dynamical instabil¬ 
ity with or without a uniform magnetic field. 


II. MODEL AND METHOD 

The Hamiltonian of the 5 = 1 BHM is given as J33ll 
T-L = —t — F 

(i,j) 7 * 

(^ 2 - 2fii ) ’ (1) 

i i 

where t is the hopping amplitude of bosons, (i,j) in the sum¬ 
mation denotes the pairs of nearest neighbors, fi is the chem¬ 
ical potential, Uq(> 0) is the on-site spin-independent repul¬ 
sion, and U 2 is the on-site spin-dependent interaction. In cold 
atom systems, the U 2 value depends on the s-wave scatter¬ 
ing length, which is specific to atom species; for example, 
U 2 > 0 (< 0) for Na (Rb) atoms. a, r , is the annihilation 
operator of a boson at site i with a spin state 7 (= 0 , ± 1 ), 
the local particle number operator n,; = a\ and 

Si = 'y2' ! y 7 ^ 7 , 7 '®*, 7 ' is the spin operator at site i where 
5 7: y corresponds to the spin-1 matrices. The square of the 
local spin operator Sf is represented in a more convenient for¬ 
mula: Sf = (5^-5^+ + Si^ + Si t -)/2 + Sf z where the ladder 
operators are defined by 5^_ = V^(a\ _ 1 di j o + 0 di,i) and 

&,+ = ., correspondingly. This formula can also be writ¬ 

ten in terms of creation and annihilation operators. 

Si — (Ay 1 1 ) “f Tli -f- Tlifi T T ^ttyofiy—1 

+2at 1 at _ 1 (Ai,o) 2 + ^(ct\fi) 2 &i,icn,-i- (2) 

The last two terms in Eq. © induce spin mixing between the 
S z = ±1 and S z = 0 states, which enriches the physics of this 
model compared with spinless models or multi-component 
models without any mixing of components. 

We first investigate the quantum dynamics of this model 
within the dynamical Gutzwiller scheme HI- The varia¬ 
tional wave function for the 5 = 1 BHM can be written as the 


direct product of superposition states at each lattice site 

'y ' fi (ttyi 5 i)|rtyi, i) > 

o 

n i,±l 

(3) 

where n^o, denotes the local Fock state deter¬ 

mined by the local number of atoms for each spin com¬ 
ponent at site i. Here the Gutzwiller parameters are nor¬ 
malized as E 7 ,n ii 7 \fi( n i,h n i,(h rii-i)\ 2 = 1. Minimizing 

(T'g|— 'HI'E'g) on the basis of the time-dependent vari¬ 
ational principle, we derive equations of motion with respect 
to these Gutzwiller parameters (47j| . The equations are ex¬ 
plicitly shown in the Appendix. Note that p(t ) in Eq. (IA2I) 
corresponds to the relative momentum between a condensate 
and a lattice for a condensate on a moving lattice or a mov¬ 
ing condensate on a stationary lattice. We introduce p{t) as 
the phase difference between particles at adjacent sites us¬ 
ing the transformation: ay 7 i —> ay 7 e JJJ( -*A (note that, t rep¬ 
resents time here). In the time-evolution calculations, we as¬ 
sume p = 0 at the initial time and the system stays in the 
ground state initially for given Uq and U 2 - The momentum 
is then increased linearly with time at the acceleration rate a: 
p(t) = at. We perform this procedure almost adiabatically by 
choosing a very small rate a = 0.005. Since loss of atoms is 
neglected in our study, the total number of particles should be 
conserved during time evolution. We ensure the number con¬ 
servation from the fact that the filling n = n\ + tiq + n-\ 
(i.e., the average particle number per site) is kept constant 
within the numerical precision. The calculated system is a 
two-dimensional lattice with a unit size L = 40 x 2 with pe¬ 
riodic boundary conditions, and we set the hopping amplitude 
t = 1 as a unit of energy. In our calculation, the sum of the 
wavefunction ([3]) is limited to a finite number of states to re¬ 
duce the number of computational tasks. We confirmed that 
the truncation does not produce any noticeable differences in 
the numerical results. 

In our calculations, velocity of a superfluid flow becomes 
quantized owing to periodic boundary conditions. However, 
the decay of a flow from the initial state to the lower wind¬ 
ing number states does not occur even in a ring geometry be¬ 
cause of the conditions we assume here, i.e., at zero tempera¬ 
ture without any thermal fluctuations. Our system therefore 
essentially becomes equivalent to the non-periodic systems 
that are generally realized in the optical lattice experiments 
to observe the dynamical instability. Actually, Mun et al. id 
reported that the observed dynamical phase-diagram can be 
quantitatively explained by the theory based on the dynamical 
Gutzwiller approximation which was developed by E. Altman 
et al. H assuming the periodic boundary condition. We also 
note that energetic instabilities like the Landau instability do 
not occur in our calculation because we keep the total energy 
in the system constant during time evolution. Therefore, there 
is no dissipation, like the phonon in the Landau instability, dis¬ 
charging energy to a heat bath such as external environment or 
a thermal component in a system. This situation is consistent 
with the experiments that observed the dynamical instability. 
As shown in Ref. d , energetic instabilities hardly appear in 


i*g>= n 

i 
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the experiments because the time scale of energetic instability 
is sufficiently longer than that of the dynamical instability at 
low temperatures where a thermal component is highly sup¬ 
pressed. 

In this paper, we discuss the instability of a superfluid flow 
by introducing two characteristic momenta: a critical momen¬ 
tum p c and a decay momentum pd. p c corresponds to the crit¬ 
ical momentum at which a superfluid flow starts to decay un¬ 
der the condition that its momentum is increased adiabatically 
from zero in an optical lattice. On the other hand, pd is the 
similar critical momentum when the momentum of a super¬ 
fluid is increased at a certain acceleration rate a as in the real 
experiment situations. A superfluid flow actually starts to de¬ 
cay drastically at pd owing to the dynamical instability during 
time evolution governed by Eq. (1 A1 b based on the dynamical 
Gutzwiller method mentioned above. From these definitions, 
Pd agrees with p c in the limit of a = 0. One can evaluate p c by 
the extrapolation using the values of pd at several acceleration 
rates, while we adopted the alternative approach mentioned 
below. 

Here we briefly explain the way we calculate p c for the 
spinless case as a simple example (ESI]- The critical mo¬ 
mentum p c is determined from the (non-dimensional) group 
velocity v(p) = p(p) sin(p) where p(p) is the density of a 
steady superfluid flowing with momentum p. The periodicity 
of v(p) reflects the structure of the lowest Bloch band in an 
optical lattice. In the framework of the Gutzwiller approxima¬ 
tion, the density p(j>) is equivalent to the condensate fraction 
Uk—p = {d' k=p dk= P ) = | (afc= P ) | 2 defined as the population of 
the state with momentum p, where k is the quasi-momentum 
of a condensate in an optical lattice. The condensate frac¬ 
tion 7ifc =p (oc t'/U q) is a monotonically decreasing function 
of p according to the effective hopping amplitude t' given by 
t' = t{d + cos p — 1 )/d where d corresponds to the dimen¬ 
sion of the system. Consequently, the group velocity v(j>) has 
a maximum at a certain momentum p = p c {< 7t/2) as p is 
increased. Beyond this p c , the effective mass, which is the 
inverse of the hopping amplitude in the tight binding model, 
becomes negative and then the sound velocity for the BHM 

becomes complex due to the formula c s = j- \J based on 
the Gutzwiller approximation [4§j], where p is the superfluid 
density and n = is the compressibility. A superfluid flow 
is unstable above p c on such a mathematical background. Fi¬ 
nally, the critical momentum p c is obtained self-consistently 
under the condition that the group velocity achieves its maxi¬ 
mum value: 


p c = arctan 


/ 

V 


n k=p=p c 




(4) 


Note that this equation is also applicable to the 5 = 1 BHM. 
We determine the phase boundary of the dynamical instability 
using this p c to remove the influence of momentum accelera¬ 
tion rate a, while the previous work employed pd 0 ]. 


HI. RESULTS AND DISCUSSION 
A. Dynamics of superfluid flow 


Figure |T| shows the time-evolution of condensate fractions 
Uk-p for the spin-1 (5 = 1 ) and spinless (5 = 0) BHM 
with filling n = 1 and on-site repulsion Uo = 10 , as func¬ 
tions of increasing momentum in time where p(t) = 0.005 1. 
In the 5 = 1 BHM, we set U 2 /U 0 = —0.3 (ferromagnetic) 
and U 2 /U 0 = 0.3 (antiferromagnetic). We choose the ground 
state as the initial state at time t = 0 for each Uq and LV- 
The total S z (= ,5V z ) in the system is conserved during 
the time evolution. From Fig. [7] nk= P gradually decreases as 
the momentum pit) is increased, and then suddenly decays 
owing to the dynamical instability. The decay momenta are 
correspondingly pd = 0.44 for U 2 /U 0 = —0.3 and 0.45 for 
U2/U0 = 0.3 in the 5 = 1 BHM, and 0.38 in the 5 = 0 BHM. 
We find that the 5 = 1 condensate persists to a larger p(t) than 
the 5 = 0 condensate, indicating some influences of the spin- 
dependent interaction included in the Hamiltonian Eq. ([]}. In¬ 
terestingly, the initial condensate fraction at p( 0 ) = 0 in the 
5=1 BHM with U 2 /U 0 = 0.3 is almost the same as that in 
the 5 = 0 model. Moreover, even in the 5=1 BHM, pd for 
U 2 /U 0 = 0.3 is slightly larger than that for U 2 /U 0 = —0.3, 
while the condensate fraction around p = pd for U 2 /U 0 = 0.3 
is apparently smaller than that for U2/U0 = —0.3. These 
results suggest that the amplitude of the condensate fraction 
does not solely determine pd, which is very consistent with 
the fact that the derivative dr ^ p p is also included in Eq. ((4} for 
determining p c . 

As we briefly mentioned in the previous section, the decay 
momentum pd inevitably becomes larger than the critical mo¬ 
mentum p c of Eq. (0} when the system parameters are equal, 
i.e., the same interaction strength, filling, and lattice geometry. 
Our previous work in Ref. [48] showed that a superfluid can 
flow stably beyond the critical momentum p c until the unsta¬ 
ble mode that causes dynamical instability fully grows. This 
retardation effect always exists as long as a finite acceleration 
of a condensate exists in calculations or experiments. We con¬ 
firmed in Fig. [j] that ^ approaches p c in the both BHMs using 
a smaller coefficient a (< 0.005) forp(f) = at. 

Next we discuss the role played by the spin mixing pro¬ 
cesses during the time evolution in the 5 = 1 BHM, which 
is governed by the third and the fourth terms on the right- 
hand side of Eq. (1A1I) in Appendix A. We focus on the an¬ 
tiferromagnetic case with U 2 /U 0 = 0.3, in which the spin 
degrees of freedom are unfrozen. In our calculations, all par¬ 
ticles are in the S z = 1 state and spins are completely frozen 
in the ferromagnetic case of U 2 < 0. Figure [2] shows the 
time-evolution of the condensate fraction rik= P and the popu¬ 
lation of each spin component n 7 /n (7 = 0 , ± 1 ) for two in¬ 
teraction strengths: (a) Uo/Uo c = 0.2 and (b) Uq/Uq c = 0 . 8 . 
Here Uq c denotes the critical interaction strength at the Mott- 
insulator transition point in the 5 = 1 BHM and Uq c = 37.9 
for U 2 /U 0 = 0.3. Note that in Fig.[2]both ri\ and n_ 1 are 
always equal owing to the initial state we choose and the con¬ 
servation of total S z . For Uo/Uq c = 0.2 shown in Fig.[2](a), 
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FIG. 1. Time evolution of condensate fraction nk= P corresponding 
to three different cases: the S = 1 Bose-Hubbard model (BHM) 
with U2/U0 = 0.3 (solid line) and U2/U0 = —0.3 (dashed line), 
and the 5 = 0 BHM (dotted line). The momentum p is increased 
almost adiabatically in proportion to time t such that p(t) = 0.005t. 
We employ the system parameters Uo = 10 and n = 1. The decay 
momenta are correspondingly pd = 0.45 for U2/U0 = 0.3 and 0.44 
for U2/U0 = -0.3 in the S' = 1 BHM, and 0.38 in the S = 0 BHM. 


the populations of the S z = ±1 states gradually decrease 
and that of the S z = 0 state increases with increasing mo¬ 
mentum, and finally all the spin components mix chaotically, 
which is accompanied by the decay of the superfluid flow. We 
also find the similar chaotic mixing of the spin components 
for Uo/Uqc = 0.8 in Fig. 0(b) after the populations of the 
S z = ±1 states have slightly increased and that of the S z = 0 
state decreased. However, the variation in spin populations 
is very small during the time evolution, suggesting that the 
spins are almost frozen in this case. We can naturally under¬ 
stand these results by noting that the third and fourth terms 
in Eq. (I A1 b make a greater contribution to spin-mixing in the 
region where Uo is sufficiently small and the amplitude of the 
Gutzwiller parameters > 2 )| 2 becomes larger. 


B. Phase diagram at unit filling 

In this subsection, we discuss how the spin mixing pro¬ 
cesses in the 5 = 1 BHM affect the critical momentum p c 
by focusing on the simple case of unit filling (i.e., n = 1 ). 
Figure[3](a) shows the dynamical phase diagram of the 5=1 
BHM with ferromagnetic interaction U 2 /U 0 = —0.3 along 
with the results of the 5 = 0 BHM. Each line represents 
the critical momentum p c as a function of interaction strength 
Uq and corresponds to the phase boundary that separates the 
stable and unstable phases. The dynamical instability occurs 
in the upper unstable region. Note that these phase bound¬ 
aries are determined via the maximum of group velocity as 
is explained in Sec ,|TT] Figure0(a) shows that the critical mo¬ 
mentum of the dynamical instability changes smoothly from 
p c = 7t/2 at Uo = 0 to p c = 0 at Uo = Uo c (i.e., the interac¬ 
tion strength at the MI transition point in the thermal equilib¬ 
rium). The critical interactions are Uq c = 33.3 for the 5=1 
BHM and U 0c = 23.3 for the 5 = 0 BHM. The cross at 


(a) (b) 



p(t)/n p(t)h t 

FIG. 2. Time evolution of condensate fraction nk= P and population 
of each spin component n 7 /n (7 = 0, ±1) in the 5=1 BHM with 
antiferromagnetic interaction U2/U0 = 0.3: (a) Uo/Uoc = 0.2 and 
(b) Uo/Uoc = 0.8. Here Uo c (= 37.9) is the repulsive interaction 
strength at the Mott-insulator transition point. We set the filling at 
n = 1. 

about p/n = 0.44 in Fig. [3] (a) on the dashed vertical line at 
Uq = 10 represents the decay momentum of a spin -1 super¬ 
fluid flow, pd, seen in Fig.[T] The apparent discrepancy be¬ 
tween this point and the phase boundary is due to the retarda¬ 
tion effect in the dynamical instability as explained in relation 
to Fig-IH 

We examine this dynamical phase diagram in more detail. 
In our calculations, all spin-1 particles with ferromagnetic in¬ 
teraction stay in the S z = 1 state, which makes the situation 
relatively simple. Therefore the spin dependent U 2 term in the 
Hamiltonian Eq. (|T|) becomes 

y^n,(ni-i). (5) 

i 

Since this form is equal to the spinless Uo term in the Hamil¬ 
tonian, the U 2 term gives just the shift in the Uq value (i.e., 
Uq —> Uq + U 2 ). In the present case of U 2 = — 0.3Uq, Uq is 
effectively reduced to 0.7?7o. As is shown in Fig. 0(b), both 
phase boundaries overlap completely when Uo is normalized 
with each Uq c - Thus the dynamical instability in the ferro¬ 
magnetic case is essentially equivalent to that in the spinless 
case. 

Next, we discuss the antiferromagnetic case. Figure[4](a) 
shows the dynamical phase diagram of the 5 = 1 BHM with 
antiferromagnetic interaction U 2 /U 0 = 0.3 along with the re¬ 
sults of the 5 = 0 BHM. We find that the two phase bound¬ 
aries are very close together for Uo < 5 (the 5=1 boundary 
is slightly below), and gradually diverge for Uq 7 ) 5. This 
divergence of the phase boundaries for Uq > 5 basically orig¬ 
inates from the difference in the Mott-transition points be- 
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FIG. 3. Dynamical phase diagrams of the S = 1 BHM with fer¬ 
romagnetic interaction U 2 /U 0 = —0.3 and the 5=0 BHM as a 
function of (a) Uo and (b) Uo/Uoc- A superfluid becomes dynam¬ 
ically unstable in the region above the phase boundary. The arrow 
in (a) corresponds to the horizontal axis in Fig.Q] and the cross in¬ 
dicates the decay momentum of the dashed line in Fig.Q] The phase 
boundaries of both models in (b) are completely identical. We set the 
filling at n = 1 . 


tween the 5 = 0 and 1 BHMs. In the strongly correlated 
regime, the probability of double occupation rij = 2 at each 
site in the 5 = 1 model is much larger than that in the 5 = 0 
model for the same interaction strength Uq owing to the for¬ 
mation of the local spin singlet state \rij,Sj : Sj, z ) = | 2 , 0, 0) 
[this formulation is defined by formula (23) in reference j38111, 
which has the energy gain — 2(72 in the (7 2 term. This en¬ 
hanced number fluctuation leads to a larger critical interaction 
strength at the Mott-transition point Uq c in the 5=1 BHM 
with unit filling. Correspondingly, the stable area of the 5=1 
model in the phase diagram grows compared with that of the 
5 = 0 system. 

In Fig. 0(b), we present the same phase diagram as a func¬ 
tion of normalized interaction strength Uq/Uq c . The phase 
boundaries of the two models are very close together for 
Uo/Uo c > 0.6, gradually diverge for Uq/Uq c 1$ 0.6, and fi¬ 
nally reach p c = n/2 at Uo/Uq c = 0. In contrast to the ferro¬ 
magnetic case shown in Fig. 0(b), the different phase bound¬ 
aries for Uq/Uoc 0-6 clearly reflect that the spin mixing 
processes included in the 5 = 1 BHM [i.e., the last two terms 
on the right hand side of Eq. (0)] play a role in the antifer¬ 
romagnetic case and influence the critical momentum of the 
dynamical instability. We examine this effect by dividing the 


phase diagram into two regions: a strongly correlated regime 
in which the phase boundaries are close together (region 1 for 
Uo/Uo c > 0 . 6 ) and another regime in which the boundaries 
diverge (region 2 for Uo/Uq c < 0 . 6 ). 

We begin by explaining region 1 in which the phase bound¬ 
aries overlap. For this purpose, we assume a system where the 
maximum number of particles per site is n max = 2 because 
the number fluctuations are greatly suppressed in a strongly 
correlated regime and the probability of rij > 3 states is neg¬ 
ligible. By noting that the U 2 term vanishes for the n ? = 1 
state, only the states, | rij, Sj, Sj tZ ) = |2,0,0) and | 2 , 2 , 77 ) 

(77 = 0 , ± 1 , ± 2 ), have non-zero energies corresponding to the 
U 2 term: -2(7 2 for | 2 , 0 , 0 ) and U 2 for \2, 2 , 77 ). The | 2 , 2 , 77 ) 
states are degenerate under the current condition without a 
magnetic field. For simplicity, we define the local spin states 
as |Sj = 0) = |2,0,0) and 15,- = 2) = ^ £ f; \2,2, V ). The 
population of the local singlet state | Sj = 0) included in the 
rij = 2 state is evaluated via the Gutzwiller parameters: 

p = _ (* G \S j =0)(S J ;=0|fr G ) _ 

0 OM Sj = 0 )(Sj = 0|tf G > + <* G |Sj = 2 ){Sj = 2|VF G ) ’ 

1(^g|5j = 0)| 2 

|(^'g|5j = 0 )| 2 + |<H/ G |5j = 2 )| 2 ’ 

§ 1 /( 0 ; 2 , 0 ) — a / 2 /( 1 , 0 ; 1)| 2 

Sn i! i+n ij0 +n i| _i=2 l/( n i,li n i,0; ?H,-l)| 2 

Furthermore, the population of the | Sj = 2) state is given by 
P 2 = 1 - P 0 . In Fig.|5](a), we show Po at both p = 0 (solid 
line) and p = p c (dashed line) as a function of Uq. Note that, 
in region 1, there is hardly any change in Po or P 2 irrespec¬ 
tive of the interaction strength and the superfluid momentum, 
reflecting the fact that the spin state becomes stationary in this 
region. This result is consistent with the slight spin variation 
seen in Fig. 0(b). The spin mixing process does not occur in 
region 1 , and consequently the spin dependent U 2 term causes 
only the shift in Uq as in the ferromagnetic case. We can thus 
understand that the phase boundaries of both the 5 = 0 and 1 
BHMs become identical when Uq is normalized by the corre¬ 
sponding Uqc as in Fig. 0(b). 

On the other hand, in region 2, the spin configurations be¬ 
come complex because the population of the rij > 3 states 
increases, and the spin mixing in the U 2 term plays a role 
(Fig. 0(b)). Here we examine how the spin degrees of free¬ 
dom influence the value of the critical momentum p c and dis¬ 
cuss the origin of the divergence of the phase boundaries in 
region 2 seen in Fig. 0(b). First, it follows from Eq. (0) that 

p c is monotonically proportional to nk= p — Pc / 

because dn ^ v < 0 for a stable superfluid flow and nfc =p > 0 . 
Furthermore, in a weakly correlated regime, the condensate 
fractions in both the 5 = 0 and 1 systems are sufficiently large 
and equally close to 1. The difference between the p c values 
of these two systems is therefore determined largely by the 
| drik=p/dp\ factor in the denominator of the above relation. 

Next we explain the influence of spins on this factor of 
\drik=p/dp\. It is generally known that the effective hop¬ 
ping amplitude of a condensate carrying the momentum p be- 
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FIG. 4. Dynamical phase diagrams of the S = 1 BHM with anti¬ 
ferromagnetic interaction U2/U0 = 0.3 and the 5 = 0 BHM as a 
function of (a) Uo and (b) Uo/Uoc. In (b), we divide the phase dia¬ 
gram into two regions: region 1 for Uo/Uoc > 0.6 and region 2 for 
Uo/Uoc < 0.6. We set the filling at n = 1. 
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FIG. 5. (a): Populations of the local spin singlet state \Sj = 0} at 
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rij > 3 states within our truncated Fock space at p = 0 (dashed line) 
and p = p c (solid line). We employ U 2 /U 0 = 0.3 and n = 1 as in 

Fig-0 


comes t' = t cos(p). The increment in momentum hence di¬ 
minishes the condensate fraction rik= P as is shown in Fig. [I] 
which simultaneously reduces the number fluctuations of a 
condensate. In the 5 = 1 BHM, this effect becomes more 
prominent thanks to the U 2 term in the Hamiltonian. Fig- 
ure[5] shows that the population of the 5) = 0) state in¬ 
creases with increasing momentum in the weakly correlated 
regime, while the population of m > 3 states decreases. This 
suggests that the number fluctuations in the 5=1 system 
are further suppressed in order to gain an energy of —2(72 
in the U 2 term. Therefore, in the weakly correlated regime, 
\dnk- p /dp\ when 5 = 1 is generally larger than that when 
5 = 0. We can confirm this fact numerically from our present 
result: 5rik=p = nk = P =0 — tifc= P = Pc is 0.045 for 5=1 with 
U 2 /U 0 = 0.3, which is nearly two times higher than 0.025 
when 5 = 0, for Uq/U c = 0.1 in both cases. Returning 
to Fig.[4](a), the difference between the phase boundaries of 
5=1 and 0 is very small in the weakly correlated regime 
owing to the small U 2 values there. However, the normalized 
phase diagram in Fig.[4](b) successfully extracts the existence 
of the spin effect on the dynamical instability of a superfluid 
flow. 


C. Phase diagrams at other fillings 

1. Commensurate case 

It is generally known that the SF-MI transition in the BHM 
strongly depends on fillings (i.e., the average number of par¬ 
ticles per site). Specifically, in the 5 = 1 BHM with antifer¬ 
romagnetic interactions, the critical interaction strength at the 
transition Uo c shows a clear dependence on the parity of fill¬ 
ings: Uqc at odd fillings is larger than that in the 5 = 0 BHM 
system, while it becomes smaller at even fillings flblfT/ll . This 
property is easily understood from the fact that the formation 
of the local singlet state to gain an energy of — 2 U 2 in the 
U 2 term in the Hamiltonian Eq. <[TJ enhances (suppresses) the 
density fluctuations at odd (even) fillings. Here we discuss 
how the parity affects the dynamical instability in the 5=1 
BHM. 

In Fig.[6] (a)-(c), the dynamical phase diagrams of the 5 = 
1 BHM for U 2 /Uii = 0.3 are given for the several different 
fillings (i.e., n =2, 3, and 4) along with the results of the 
5 = 0 model. From these figures and Fig. [4] (a), we find that 
the influence of the parity clearly appears in the dynamical 
phase diagrams. The stable areas of the 5 = 1 model basically 
grow (shrink) at even (odd) fillings compared with the 5 = 0 
model, which reflects the corresponding increase (decrease) 
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FIG. 6 . Dynamical phase diagrams in the S = 1 and 5 = 0 
BHMs for three different fillings: (a) n = 2, (b) n = 3, and (c) 
n = 4. In (d), the results of (a) are presented as a function of Uo/Uo c . 
U2/U0 is fixed at 0.3 in all cases for the S' = 1 model. Critical 
interaction strengths are Uoc = 39.6 ( n = 2), 55.2 (n = 3), and 
70.7 (n = 4) for 5 = 0, and U 0c = 11.2 (n = 2), 68.45 (n = 3), 
and 18.0 (n = 4) for 5 = 1. 


of Uq c . With n = 3 shown in Fig.[ 6 ](b), however, the stable 
area clearly decreases in the weakly correlated regime. As 
we pointed out for unit filling in the previous subsection, the 
unfrozen spins that prefer to form the local singlet state greatly 
suppress the density fluctuations and make the superfluid flow 
unstable in the weakly correlated regime. We have confirmed 
this effect more clearly for n = 3 filling. This result indicates 
that the spin mixing process has a greater influence at larger 
fillings. 

In Fig.[ 6 ](d), we provide a dynamical phase diagram for 
n = 2 filling as a function of the normalized interaction 
Uo/Uoc The 5=1 phase boundary is located above the 
5 = 0 curve over the entire interaction range, which is in con¬ 
trast to the phase diagram for unit filling shown in Fig.[4](b). 
This suggests that the unfrozen spins, which prefer the local 
singlet states, stabilize the superfluid flow. We have also found 
this tendency with n = 4 filling as a characteristic of the dy¬ 
namical instability at even fillings. 


2. Incommensurate case 

In a system with an incommensurate filling, the SF-MI tran¬ 
sition does not occur because of the extra particles deviating 
from the commensurate filling. Polkovnikov et al. calcu¬ 
lated the dynamical phase diagram of the 5 = 0 BHM with 
incommensurate fillings based on the Gutzwiller approxima¬ 
tion in Ref. Ht]. They clarified that the critical momentum 
p c has a minimum value at a certain interaction strength and 
then asymptotically approaches 7r/2 with increase in interac¬ 
tion strength. The superfluidity of the extra particles becomes 
highly robust in the strongly interacting regime, i.e., the su¬ 
perfluidity recovers owing to the repulsive interaction. The 
minimum of p c in the dynamical phase diagram therefore rep- 


FIG. 7. Dynamical phase diagram for incommensurate fillings: (a) 
n = 0.8, (b) n = 1.2, (c) n = 1.8, and (d) n = 2.2. We employ 
U 2 /U 0 = 0.3 in all cases for the 5 = 1 model. 


resents the crossover between weakly and strongly interacting 
regimes. Here we analyze this tendency in the 5 = 1 BHM 
with the antiferromagnetic interaction U 2 /U 0 = 0.3. 

Figure|7] shows the dynamical phase diagrams of the 5=1 
BHM with the filling factors deviating slightly from n = 1 
and 2, along with diagrams of the 5 = 0 BHM. From these 
figures, the dynamical phase diagrams of the 5=1 BHM at 
incommensurate fillings agree qualitatively with those of the 
spinless 5 = 0 model. However, we still find the influence of 
the parity, which we have seen in the phase boundaries for the 
commensurate cases. In Fig. 0(a) and (b) where the fillings 
are close to n = 1 , the critical momentum reaches its mini¬ 
mum value at a larger interaction strength in the 5 = 1 BHM 
in comparison with the 5 = 0 results. On the other hand, in 
Fig.[7](c) and (d) where the fillings are close to n = 2, the 
minimum point for 5 = 1 is apparently smaller than that for 
5 = 0. This behavior can be roughly understood in terms 
of whether the formation of the local spin singlet state con¬ 
ceals or accentuates the extra particles. As mentioned in the 
commensurate case, the formation of the local singlet state 
enhances (suppresses) the density fluctuation of a condensate 
with the odd (even) fillings. The intense density fluctuation 
of a condensate conceals the effect of the extra particles on 
the left of the minimum points while the suppressed density 
fluctuation accentuates the extra particles on the right side. 
Therefore, with the fillings close to n = 1, the extra particles 
are more concealed due to the formation of the local spin sin¬ 
glet state, and the minimum point for 5 = 1 slides to the right 
compared to that of 5 = 0, while the extra particles are more 
accentuated due to the formation of the local spin singlet state 
and the minimum point for 5 = 1 slides to the left compared 
with that of 5 = 0 with the fillings close to n = 2. Fur¬ 
thermore, we see that there is particle-hole symmetry in the 
dynamical phase diagrams with incommensurate fillings by 
noting the consistency between the n = 0.8 and 1.2 results, 
and also between the n = 1.8 and 2.2 results. 
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site i site i 


FIG. 8. Density modulation associated with the dynamical insta¬ 
bility with no magnetic field, (a): Density distributions of S z = 1 
(triangles), S z = 0 (circles), and S z = — 1 (squares) components 
over the lattice sites. The density distributions of S z = ±1 are fully 
identical, (b): Deviation of densities from the mean values: particle 
number (light line) and magnetization (dark line). The parameters 
are fixed at n = 1, Uo = 10, and U 2 /U 0 = 0.3. These results corre¬ 
spond to the point atp/n = 0.46 after the decay at Pd/n = 0.45 on 
the solid line in Fig.Q] 


D. Density modulation 

Finally, we discuss the density and spin modulation asso¬ 
ciated with the dynamical instability. A spinless superfluid 
flow in an optical lattice exhibits a density modulation as a 
precursor to dynamical instability Oil)]. This is a mani¬ 
festation of unstable collective excitation modes as a seed of 
the dynamical instability, and it depends strongly on the inter¬ 
action strength or the acceleration rate of the condensate mo¬ 
mentum 114811 . This collective excitation, which involves a lot 
of physical information, is significant in terms of understand¬ 
ing the dynamical instability. Here we examine whether the 
occurrence of spin modulation is associated with the dynam¬ 
ical instability of the spin-1 condensate in the 5 = 1 BHM. 
We again assume an antiferromagnetic system with a spin- 
dependent interaction strength of U 2 /U 0 = 0.3. 

Figure[ 8 ](a) shows the density distributions of S z = 0, ±1 
components after the dynamical decay of a condensate. The 
density modulation develops sufficiently at this momentum. 
We find that the density modulations of the S z = ±1 compo¬ 
nents develop in unison, and that of the S z = 0 component 
develops independent of those modulations. This result indi¬ 


FIG. 9. Density modulation associated with dynamical instability 
under a uniform magnetic field, (a): Density distributions of S z = 1 
(triangles), S z = 0 (circles), and S z = — 1 (squares) components 
over the lattice sites, (b): Deviation of densities from the mean val¬ 
ues: particle number (light line) and magnetization (dark line). The 
parameters are n = 1, Uo = 10, and U2/U0 = 0.3 as in Fig[8] 
A uniform magnetic field is used to adjust the initial populations of 
the spin components to m : n_i ~ 7 : 3 under the condition of 
p — 0. These results correspond to the point at p/n = 0.44 after 
the decay at pd/it = 0.427 as the momentum is increased such that 
p(t) = 0.005f. 


cates that small spatial fluctuations in a condensate grow in¬ 
dependently in S z = ±1 components and the S z - 0 compo¬ 
nent. This is because the components of S z = ±1 are equiv¬ 
alent in Eq. (IA1I) within the mean-field approximation where 
the density modulations of S z = ±1 components develop in 
unison. Figure[ 8 ](b) shows the total density distribution and 
magnetization distribution at the same momentum. There is 
no spin modulation of the S z components while the total den¬ 
sity modulation develops intensely. This reflects the consis¬ 
tent development of the modulations of the S z = ±1 compo¬ 
nents. Therefore, the spin modulation occurs only within the 
xy plane. 

Next, let us examine the spin modulation in a system with 
a uniform magnetic field in the z direction. Here, the S z = 
±1 populations become imbalanced and we adjust them to 
n\ : n-i ~ 7 : 3. Figure[9](a) shows the density distribu¬ 
tions of S z = 0, ±1 components after the dynamical decay 
of a condensate. A magnetic field is applied to the system 
only at initial state p = 0, but the initial state is stable since 
the total S z in the system is conserved. A significant differ- 
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ence from non-magnetic case is that the density modulations 
of the S z = ±1 components develop independently, namely, 
small spatial fluctuations in a condensate grow independently 
in S z = 1 and S z = — 1 components. This result indicates 
that S z = ±1 components decay independently only if there 
is difference between the component populations. As a re¬ 
sult, the spin modulation of S~ components occurs as shown 
in Fig.[9](b). 


IV. SUMMARY 

In this study, we analyzed the dynamical instability of a 
superfluid flow in the 5 = 1 BHM using the Gutzwiller ap¬ 
proximations. Time evolutional calculations revealed that the 
superfluid flow of the spin -1 condensate decays at a different 
critical momentum from the spinless model when the inter¬ 
action strength is the same, which is due to spin-dependent 
interactions. Furthermore, we obtained the dynamical phase 
diagrams of both the 5 = 1 and spinless 5 = 0 BHMs and 
discussed their differences. With a ferromagnetic interaction 
U 2 < 0, the phase diagram of the spin-polarized 5=1 BHM 
becomes essentially the same as the diagram of the spinless 
BHM because we can appropriately renormalize the interac¬ 
tions. On the other hand, with an antiferromagnetic interac¬ 
tion U 2 > 0, the dynamical phase diagrams of the 5=1 
BHM differ fundamentally from the spinless model and shed 


light on the influence of the spin mixing process between the 
5 = 1 bosons. We discussed in detail the important role of the 
formation of the local singlet state in the dynamical instability 
and the SF-MI transition in the 5 = 1 BHM. Our systemati¬ 
cal study also showed that the phase diagram strongly depends 
on the average number of particles per site, in particular, the 
even-odd parity. We finally discussed the density modulation 
and the spin modulation associated with the dynamical insta¬ 
bility. We found that the anisotropy of the spin modulation de¬ 
pends on whether or not a uniform magnetic field is present. 
This suggests that the spin modulation is highly sensitive to 
the imbalance in the spin components generated by a uniform 
magnetic field. 
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Appendix A: The equations of spin-1 Gutzwiller parameters 

The equations of motion for the Gutzwiller parameters in 
the 5 = 1 BHM are given as 


n j, 0 i n j,l) — 2 n j( n 3 — 

+ - 2n j,i n j,-i + n j,-i - n jA ~ n j-i + 2n j,i n j,o + 

YC 2 \ Jfij 1 (fij 0 T 1 )(fijo -p 2 )iij _ | /'.y (Tij _ 1 1, vt.^0 T" 2, 1) 

+U 2 yJ (raj, 1 + l)ttj,o(jtj,o - l)( n i,-i + ,-1 + l,nj )0 - 2,71^1 + 1 ) 

^ (\AC. — i I ;/(. — 1 I? Y \j 1 Y 1 fj ( Tij ,—1 4- 1, 

-tz - 1, nj, i)ipj,o Y \Au,o Y 

-tz Y Y _i, n ji0) n jt i Y l)^,i) , 


(Al) 


where 


i’j, 1 — 7 E V n j+ 1. 1 + ^-fj+i( n j+i-iy n j+i,o,nj+i t i)fj+i(nj + i j -i,nj + ifi,nj + i } i Y l)e* p ^ 

n j + l,0 
n j + 1,±1 

H \J n j~ 1,1 + n j-l,0i 1,0 5 n j-l,l + l) e 


nj- 1,0 
rij- i,±i 


+tE E 7^7TT/ j+r 1 ? W'j+T, O 7 (jlj+T, — 1 5 Ch Wj+r, 1 H - !)• 


r n j+T, 0 
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4 >j, o = 7 E \/ n i+i,o + 'i-fj+i(nj+i-i,nj + i } o,nj + i t i)fj + i(nj + i j -i,nj + ifi + 1 , 'rij+i,i)e lp< '^ 

n J + l,0 
n J + l,±l 

+ 7 V n j- i,o + + l, n j-i,i) e 


-ip(t) 


Tlj — 1,0 
Tl'j — 1, i 1 


+tE E \/^'i+' 7 ',0 H“ (^j+r,— 1 ? ^j+TjO? ^.7 +r, 1) /j+r {^j +t ,'—1 ? ^-fr,0 “I” 1? 


r n j+T,o 
^i+T,±i 


V’jj-i — 7 E V n J + 1 11 -i + l/j*+i( n J+i,-i! n j+i,o,«j+i,i)/i+i(«j+i,-i + l,«j+i,o,«j+i,i) e * P< '^ 

n J + l ,0 

n J + l,±l 

+ 7 E \/ n i-L-i + -i,«i-i,o,^j-i,i)/j-i(«i-i,-i + l,nj_i, 0) nj_i,i)e _ip(t) 


n j-l ,0 
n i—i,±i 


+tE E \/^+T, —1 “1” (jlj+T, — 1 ? ^J+TjO? ^'J+T,l)/ 1 7-|-T (P'j+T,— 1 H - 1? ^J+TjO? ^i+T,l ) • 


T n i+r,0 
^lj + r,±l 


(A2) 


Here j + 1, j — 1, and j + r denote (ji +1, J 2 ), (ji — I.J 2 ), summation runs over the nearest neighbors of site j in the 
and (ji,j 2 + t) respectively, where ji is the site index of the orthogonal direction, and z is the number of adjacent sites in 
flow direction and j -2 is that of the orthogonal direction. The the lattice. 
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